Mon. Not. R. Astron. Soc. OOP. [TlfTSl C2008') Printed 8 February 2008 (MN I*TeX style file v2.2) 



Planetesimal and Gas Dynamics in Binaries 



S.-J. Paardekooper^*, P. Thebault^'^ and G. Mellema^ 

^ DAMTP, University of Cambridge, Wilberforce Road, Cambridge CBS OWA, United Kingdom 
^Stockholm Observatory, Albanova Universitetcentrum, SE-10691 Stockholm, Sweden 
QQ ' "^LESIA, Observatoire de Paris, Section de Meudon, F-92195 Meudon Principal Cedex, France 



Draft version 8 February 2008 



ABSTRACT 

Observations of extrasolar planets reveal that planets can be found in close binary sys- 
tems, where the semi-major axis of the binary orbit is less than 20 AU. The existence of 
these planets challenges planet formation theory because the strong gravitational per- 
turbations due to the companion increase encounter velocities between planetesimals 
and make it difficult for them to grow through accreting collisions. We study planetes- 
imal encounter velocities in binary systems, where the planetesimals are embedded in 
a circumprimary gas disc that is allowed to evolve under influence of the gravitational 
perturbations of the companion star. We use the RODEO method to evolve the ver- 
tically integrated Navier-Stokes equations for the gas disc. Embedded within this disc 
is a population of planetesimals of various sizes, that evolve under influence of the 
gravitational forces of both stars and friction with the gas. The equations of motion 
for the planetesimals are integrated using a 4**^ order symplectic algorithm. We find 
that the encounter velocities between planetesimals of different size strongly depend 
on the gas disc eccentricity. Depending on the amount of wave damping, we find two 
possible states of the gas disc: a quiet state, where the disc eccentricity reaches a 
steady state that is determined by the forcing of the binary, for which the encounter 
velocities do not differ more than a factor of 2 from the case of a circular gas disc, 
and an excited state, for which the gas disc obtains a large free eccentricity, which 
drives up the encounter velocities more substantially. In both cases, inclusion of the 
full gas dynamics increases the encounter velocity compared to the case of a static, 
circular gas disc. Full numerical parameter exploration is still impossible, but we de- 
rive analytical formulae to estimate encounter velocities between bodies of different 
sizes given the gas disc eccentricity. The gas dynamical evolution of a protoplanetary 
disc in a binary system tends to make planetesimal accretion even more difficult than 
in a static, axisymmetric gas disc. 
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1 INTRODUCTION 

Planetary formation in binary systems is an issue which has 
gained considerable attention in recent years, fueled by the 
discovery of more than 40 extrasolar planets in multiple sys- 
tems, making u p over 20% of the total number of exoplanets 
detected so far (jRaghavan et al. I I2OO6I : iDesidera fc Barbieril 
l2007i ). Note, however, that the vast majority of these plan- 
ets have been detected in wide binaries, with separations 
^ 100 AU, for which the binarity of the system is proba- 
bly only weakly felt at the location of the planets (most of 
which have semi-major axis Op < 3 AU). Nevertheless, a 
handful! of planets have been identified in "tight" binaries 
with separations of typically ~ 20 AU: HD 41004, Gliese 86 
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and perhaps the most famous, 7 Cephei. For these systems 
there is little doubt that the presence of the companion must 
have played a crucial role. 



The first studies of the planets-in-binaries issue were 
devoted to the problem of long term orbital stability. In 
this field, the reference work pr obably remains the analy i- 
cal and numerical exploration bv lHolman fc Wiegerd l|l999^ ■ 
which enabled these authors to derive semi-empirical for- 
mulae for the critical semi-major axis flcrit beyond which a 
planetary orbit becomes unstable in a binary of mass ra- 
tio n = m2/(mi -I- m2), separation Ob and eccentricity Cb. 
Since then other numerical studies further explored this is- 
sue, and the criteria for orbital stability in binaries are now 
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fairly well constrained (e.g . 
l2006l : lMudrvk fc Wulliooa ) 



David et al.ll2003l : iFatuzio et al.l 1.2 Planetesimal impact velocities in binaries 



The problem of the binarity's influence on planet 
formation is a different issue. In the now "standard" core- 
accretion scenario, planets are believed to appear through 
a succession of different stages, leading in a complex and 
still only partially understood way from small grains con- 
densed in the initial neb ula to fully formed planetary ob- 
jects (e.g iLissaueJ Il993l ). The way each of these stages 
can be affected by the presence of a companion star re- 
quires in principle a full study in itself. In this respect, two 
stages have been investigated in particular. One of them is 
the final mut ual assembly of protoplanetary embryos into 



toplan( 

planets (e.g. Chambers et al. 20021 : iQuintana et al.l l2007l : 



iHaghighipour fc Ravmondll2007l l. This phase is almost en- 
tirely governed by gravitation, more exactly by the cou- 
pling between mutual gravitational pertubartions among 
the embryos and the external pull of the companion star. 
lOuintana et al.l (|2007l l showed that, as a first approxima- 
tion, the final stages can proceed unimpeded at distance ap 
if the binary's periastron ah exceeds ~ 5ap. 



1.1 Planetesimal accretion 

The other stage that has been extensively studied, and the 
one which we focus on in the present study, is the earlier 
phase leading from kilometre-sized planetesimals to 100- 
500 km sized embryos. This planetesimal accretion phase 
should in principle be much more sensitive to the perturb- 
ing presence of the companion star. Indeed, "normal" plan- 
etesimal accumulation around a single star should proceed 
in a dynamically cold environment, where average mutual 
encounter velocities (Av) are of the order of planetesimal 
surface escape velocities Vcsc, i.e. only a few ms~^ for km- 
sized bodies. In such a quiet swarm, a few isolated seed ob- 
jects, initially just a little bigger than the rest of the swarm, 
can grow very quickly, through the so-called runaway growt h 
mode (e.g lGreenberg et al.lll978l : IWetheriU fc Stewart! 19891 ). 
This very fast and efficient mechanism takes place because 
these bodies' coUisional cross section ctcoii is amplified, with 
respect to their geometrical cross section (Jgeom, by a gravi- 
tational focusing factor: 



f7"coll 
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1 + 



WcscR 



(1) 



where "R" refers to the seed object and "r" to the bulk of the 
planetesimal population. This growth is self-accelerating: 
faster growth leads to higher amplification factor, which in 
its turn leads to faster growth, and so on. As can be seen 
from Eq. ([T}, runaway growth is very sensitive to (Av). In a 
binary system, this growth mode is obviously at risk: if the 
companion is efficient enough in stirring up the planetesimal 
swarm, up to velocities exceeding VeBc{R), then runaway ac- 
cretion could easily be stopped. Further stirring could even 
lead to total stop of any form of accretion, by increasing 
(Av) above the threshold velocity Uoro, for which all impacts 
lead to mass erosion instead of accretion. 



^ In tiiis respect, it is worth noticing that all planets discovered 
so far in binary systems are reassuringly on stable orbits 



This possible effect of a stellar companion on the impact 
velocity distribution was explored in several previous works. 
The first studies dealt with the simplified case where only 
the gr avitational potential of the two stars is tak en into ac- 
count iHeppenheime i<1978l : IWhitmire et al.ll 19981 ). However, 
within current planet formation scenarios, it is very likely 
that planetesimal accretion takes place while vast quantities 
of primordial gas are left in the system. The complex coupled 
effects between binary secular perturbations and gaseous 
friction on the planetes imals were numerical l y exp lored by 
iMarzari fc Scholll l|2000l ) and iThebault et"al] ^20041 ) for the 
specific cases of the g-Centauri and 7 Cephei systems, and 
more recently by IThebault et ah I l|2006l ) for the general case 
of any binary configuration. The latter study in particular 
showed that, while gas drag is able, through strong orbital 
alignement, to counteract the velocity stirring effect of an 
external perturber for impacts between equal-sized bodies, 
it tends to increase (Av) between objects of different sizes. 
This is because gaseous friction effects are size dependent. 
The main conclusion was that even relatively wide binaries 
with 10 ^ Ob ^ 50 AU could have a strong accretion inhibit- 
ing effect on a sw arm of colliding kilom etre-sized planetesi- 
mals (see Fig.8 of IThebault et al.ll2006l ). 

However, because of numerical contraints, these stud- 
ies were carried out assuming a perfectly axisymetric gas 
disc. This simple approximation is unlikely to be realistic. 
Indeed, the gas should also "feel" the companion star's per- 
turbations and depart from an axisymmetric structure with 
circular streamlines. To what extent does the gas disc re- 
spond to the secular perturbations? Will gas streamlines 
align with planetesimal forced orbits, thus strongly reducing 
gas drag effects, or will they significantly depart from forced 
orbits, thus playing a major role for planetesimal dynam- 
ics? More realistic models of the planetesimal-|-gas evolu- 
tion are definitely needed. A pioneering attempt at mod- 
eling coup led planetesima l /gas systems was recently per- 
formed by ICiecielag et al.l l|2007h . This study was however 
restricted to a circular binary, which is the most simple 
to numerically tackle, because in this case the spiral wave 
pattern that is excited in the gas disc is stationary, which 
removes the need for following the evolution of the gas 
numerically. Note, however, that possible gas eccentricity 
;rowth throu g h the effects of the 3:1 Li ndblad resonance 
Lubowl I1991I : iGoodchild fc Ogilvid I2OO6D is not included 
using this method, nor are any eff ects of the viscous over- 
stabhlity (K ato 1978; L atter fc Og ilvic 2006). Furthermore, 
ICiecielag et al.l (|2007l ) focused on following semi-major axis 
and eccentricity evolutions, without computing encounter 
velocities estimates, which is, as previously discussed, the 
crucial parameter for investigating the accretional evolu- 
tion of the swarm. Kley fc Nelson (2007.) performed a full 
scale planetesimal-|-gas study in eccentric binaries, includ- 
ing the gravitational pull of the gas disc on planetesimals. 
However, the complexity of their numerical code made that 
only a limited number of planetesimals could be followed 
(less than a hundred), which is not enough to derive relative 
velocity statistics. Basically, including the gravitational pull 
of the gas for ~ 10000 particles is equivalent to including 
self-gravity of the gas, which is known to be very expensive 
numerically. See Sect. |6] for a discussion. 
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1.3 Paper outline 

We aim to extend tiie study of Ixiiebault et alj (|2006l ) by 
including tlie efltects of an evolving gas disc. Solving the 
hydrodynamical equations of motion for the gas is compu- 
tationally expensive compared to solving the equations of 
motion for the planetesimals, which makes it impossible to 
perfor m a parameter study as presented in iThebault et alJ 
l|2006t ). Instead, we focus on two binary configurations only, 
one corresponding to a tight binary of separation ab =10 AU 
with Ch = 0.3, and the other to the specific case of 7 Cephei, 
and point out the general trends induced by dynamical evo- 
lution of the gas. 

The plan of this paper is as follows. In Sect.[2]we discuss 
our model and the numerical method. In Sect. [3] we present 
some simple analytical considerations on the problem, after 
which in Sect . [4] we test our method by reproducing some of 
the resuhs of lfhebault et all (|2006h . We come to the results 
in Sect.O which we subsequently discuss in Sect.|Sl and we 
sum up our conclusions in Sect. [7] 



2 MODEL 

Our model consists of three distinct components: the binary 
system, the planetesimal disc and the gas disc. All orbits are 
assumed to be coplanar. 



2.1 Particles 

The orbits of the binary stars of masses Mp and Mb (for the 
primary and the secondary, respectively) and the planetesi- 
mals are integrated using a fourth order symplectic method. 
The binary stars feel only the gravitational force from each 
other, while the planetesimals in addition feel friction with 
the gas (see Eq. (ITS)) '): 



fdr 



-K\Av\ Av, 



(2) 



where Av denotes the velocity difference between the plan- 
etesimals and the gas, and the constant of proportionality 
K is given by: 

3pgCd 



K 



SppiR ' 



(3) 



in which pg is the gas density, ppi is the planetesimal in- 
ternal density and R is the radius of the planetesimal. We 
consider spherical bodies only, whic h make s the drag coef- 
ficient Cd ^ 0.4. Following .ThcbauH et al.l (|2006l ). we take 
ppi = 3 g cm~^. 



2.2 Gas disc 

The gas disc is modeled as a two-dimensional, viscous ac- 
cretion disc. We will work in a cylindrical coordinate frame 
(r, 95) with the primary star in the origin. Note that this is 
not an inertial frame. It is convenient to work with units in 
which GAIp = Ob = 1. In these units, the orbital frequency 
of the binary is unity. 

The disc extends radially from r — Q.025ab to r = 0.4ab. 
The equation of state is taken to be locally isothermal: 



where p is the (vertically integrated) pressure, S is the sur- 
face density and Cs is the sound speed, which is related to 
the disc thickness H and the Kepler frequency Qk through 



(5) 



In all simulations presented here, we have used H = 0.05r. 

We do not consider self-gravity for the gas, to keep the 
problem tractable computationally. Note, however, that for 
globally eccentric gas discs self-gravity may be an important 
effect (see Sect.[51). 

The surface density follows a power law, initially, with 
index —7/4 or —1/2. Because we do not include self- gravity, 
the density can be scaled arbitrarily, as far as the hydro- 
dynamical evolution is concerned. We normalize the surface 
density so that pg = 1.4 • lO^'"* g cm^^ at 1 AU. This is 
consistent with the Minimum Mass Solar Nebula (MMSN) 
of Havashi (1981). 

The disc is initialized to have equilibrium rotation 
speeds (slightly sub-Keplerian due to the pressure gradient), 
with zero radial velocity. The viscosity is modeled using the 
Q-prescription, with a = 0.004 , which is a standard value 
for si mulations of protoplanetary discs (see for example iKlevI 
Il999f) . 



2.3 Numerical method 



We 



drody- 
This 



P ■ 



(4) 



use the RO DEO method for solving the 
namic equations (jPaardekooper fc Mellema I2OC 
method has been used to simulate planets embed- 
ded in protoplanetary di s cs con sisting of gas and dust 
(|Paardekooper fc Mellerna 2006ah. and wa s compared to 
other methods in Ide Val-Borro et al.l l|2006l ). The inclusion 
of par ticles into the method was discussed in IPaardekooperl 
(|2007l ). 

We monitor encounters between planetesimals using the 
procedure which is standard in such studies, i.e., assigning 
an inflated radius Rind to every particle. When two inflated 
particles intersect, we measure their relative velocity Av. 
This gives us a distribution in space and time of collision 
velocities. For a planetesimal disc consisting of 10000 par- 
ticles, RinC ~ 10~^ AU resulted in good enough collision 
statistics, without introducing an artificial bias in relative 
velocity estimates exceeding ~ 10ms~^. 

We use non-re flecting boundary conditions th rough- 
out this paper fsee iPaardekooper fc Mellemall2006bl ). This 
way, we avoid spurious reflections of spiral waves off the 
grid boundaries. We vary the flux limiter, which basi- 
cally controls the amount of numerical wave damping, be- 
tween the minmod l imiter and the soft superbee lim iter 
(see Appendix A and IPaardekooper fc Mellemal (|2006bl ) for 
their deflnition). We also included a kinematic a- viscosity, 
parametrized by q = 0.004. The adopted grid size is 
{nr,n^) = (192, 448), yielding a resolution Ar = 0.002 ttb. 



3 ANALYTICAL CONSIDERATIONS 

In linear theory, the response of the gas disc to a gravita- 
tional perturbation occurs only at Lindblad and corotation 
resonances. At Lindblad resonances, the perturbations take 
the f orm of a trailing spriral waves (|Goldreich fc Tremaind 
I1979I ) with pattern speed il;,™ . For a perturber on a circular 
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orbit, all modes have a pattern speed equal to the orbital 
frequency of the perturber fib • A key result is that even for 
a very close binary of ah = 10 AU the orbital frequency in 
the planet-forming region (say 1 AU) is much larger than the 
pattern speed of the perturbation, the effect of which is that 
planetesimals at 1 AU will sweep through the spiral pattern 
very quickly, and experience many short encounters with 
the density waves. In a recent study, ICiecielag et"al] (|2007l ) 
showed that only the very small planetesimals (< 100 m) 
are significantly affected by the waves. 

Modes of low azimuthal wavenumber are potentially 
much more important in perturbing the orbits of planetesi- 
mals. In this section, we look at the effect of a global m = 1 
perturbation in the gas disc, or, in other words, an eccentric 
gas disc, on the planetesimal swarm. We do not specify the 
source of the gas eccentricity, which may be the eccentric 
binary companion, the viscous overst ability (p<atQ,1978i '). or 
an eccentric resonance (|Lubowlll99ll ). 



3.1 No-drag limit 

iGoodchild fc Qgilvid { 20061 ) derived the governing equation 
for the secular evolution of the complex eccentricity E of 
a gas disc subject to an external perturbing potential on a 
circular orbit. Two straightforward modifications allow us to 
use their result on our planetesimal swarm in the limit of no 
gas drag: we let the pressure p — > and introduce an m = 
1 component to the perturbing potential. The eccentricity 
equation then reads: 



iE_d_ 

r dr 



I 8^2 

dr 



+ 



d_ 

dr 



(6) 



where Q. is the angular velocity and $2 and ^'2 are the ax- 
isymmetric and the m = 1 component of the perturbing 
potential, respectively. 

The perturbing potential can be expanded in terms of 
a Fourier cosine series in (f>, and, because we are interested 
in secular perturbations only, time averaged. We only need 
to retain t he first two terms in the Four ier series to lowest 
order in r (jAugereau fc Papaloizoul [20041 '): 
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(7) 
(8) 



where quantities with subscript "b" refer to the binary com- 
panion. 

With above expressions for the perturbing potential, 
Eq. (O can be integrated to 
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where is the secondary to primary mass ratio and Pb is the 
orbital period of the binary, and we have taken E{r, 0) — 0. 



This expression is equivalent to 
e{r,t) = 2ef(r) 



U t 
sm I — =— 
2 Pb 



tan (tu) — 
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ef(r-) 



sin U 
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(10) 
(11) 

(12) 
(13) 



When we replace r with a this corresponds exactly to the 
result given by non-linear secular perturbation theory. The 
eccentricity oscillates around the forced eccentricity et , with 
a spatial frequency that increases with time. In the first 
stages of the system's evolution, strong orbital phasing be- 
tween neighbouring objects (see Eq. Hi}) prevents orbital 
crossing. However, at some point (which depends on the ra- 
dial location within the system) the frequency becomes high 
enough for planetesimals of high eccentricity to encounter 
planetesimals of low eccentricity, resulting in orbital cross- 
ing and mutual enco unters at very high en counter velocities 
(for more details, see iThebault et al]|2006l) . 



3.2 Linear drag 

Including a linear drag law into the equation of motion 
Av 

Tf 



(14) 



where An is the velocity difference between gas and plan- 
etesimals and Tf is the friction time scale, leads to an extra 
term in Eq. (|6]) 

-flii (£-£,), (15) 

where Eg is the eccentricity of the gas. The eccentricity evo- 
lution of the planetesimals now consists of a damped oscil- 
lation: 



E{r,t) = 



15ebgb(^) Hrf -hl6iSg (1-eg)'/' 



exp 



Hi -elf 



-3/2 



The presence of the damping term (the second term in 
square brackets) ensures that the system will evolve towards 
a stationary state, with 



2\5/2 



Eoir) 



15ebgb nn + WiEg (l - eg 

12gb(^)'(l-eg)nrf + 16i(l-eg)^/^' 



(17) 



The important thing to notice is that Eq in general depends 
on Tf, which means that particles of different sizes (leading 
to different values of rf) will evolve towards a different ec- 
centricity distribution. This differential orbital phasing will 
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Figure 1. Axisymmetric gas disc case (eg = 0). Equilibrium 
eccentricity (solid curves) and zu (dashed curves) distributions 
for et, = 0.3, a^, = 10 AU and qt, = 0.5 (tight binary case) for 
planetesimals of 1 km (black curves) and 5 km (gray curves), as 
given by Eqs. I I28I I and II29I I. Also shown is the forced eccentricity 
ef (black dotted line). 



lead to very high encounter velocities between planetesimals 
of different sizes. Note, however, that for the special case of 
Eg — ef there is no size dependence in the eccentricity evo- 
lution. In this special case, all planetesimal orbits will be 
damped towards E = et, regardless of their size. For this 
case, we expect small encounter velocities. 



3.3 Non-linear drag 

A linear drag force is not appropriate for planetesimals. For 
large bodies, a non-linear drag force is more appropriate: 



/drag = -K I A?)| Al;, 

where 

3Pg 



K 



(18) 



(19) 



where pg is the gas density, ppi is the internal density of the 
planetesimals, and R is the planetesimal radius. We have 
assumed the planetesimals to be spherical. For this drag 
law, the drag term Eq. (fTS)) gets replaced by 



-VSKrn^ \E - Eg] {E - Eg) . 



(20) 



In this case it is not feasible to obtain an expression for the 
time evolution of the eccentricity, but we can still look at 
the stationary solution, which is given by: 
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z + iIC\z\ z + Eg — ef 
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Then it is straightforward to show that in this case 



Eq = ef + {Eg - ef 
where 

Z = Jl + ^K? \Eg 



Z + \ 



■ Ef 



(23) 



(24) 



(25) 



All size-dependence is contained in Z, and it is clear that 
again for the special case Eg = ef there is no differential 
orbital phasing. Also for this case we expect small encounter 
velocities. 

For completeness, we write out the magnitude eo and 
the argument tjjq of the complex equilibrium eccentricity Eq: 



2 2e^ + (^-l)e; 
eo — 



2^2Z - 2e{I{Eg 



Z + 1 



(26) 



tan(ci7o) 



{Z - l)I(Eg) - V2Z^{ei - 7^(£g)) 
{Z - l)I{Eg) - ^/2Z -2 (ef + TZ(Eg)^ 



(27) 



where TZ{Eg) and 1{Eg) denote respectively the real and 
imaginary part of Eg and eg its absolute value. The last two 
expressions lead to two interesting observations. When the 
gas disc eccentricity vector has a significant positive imagi- 
nary part, the eccentricity of the particles is reduced. This 
happens when the gas disc is not aligned (or anti-aligned) 
with the binary orbit. For the special case Z = Z and 
Eg = iet we have Eo = 0. Though this may seem a patho- 
logical case, we will see that this behavior actually shows up 
in the numerical simulations. 

The case of a circular gas disc (defined as having eg = 
0, but not necessarily being axi-symmetric) appears to be 
particularly simple: 



eo 



ef 



tan(ti7o) 



1 + Z 

Z-\ 



(28) 
(29) 



In the inner regions of the disc, where gas drag is very strong 
and therefore Z ^ 1, the eccentricity will be much smaller 
than ef, while the longitude of periastron will make a 90- 
degree angle with that of the binary. This effect can be 
clearly seen in Fig. 6 of iThebauh et all l|2006l) (also see our 
Fig. As Z ^ 1, indicating less friction (either due to 
a lower gas density or a larger particle size), particles tend 
to adopt the forced eccentricity and align with the binary. 
Note, however, that this analysis does not take into account 
the oscillatory behavior of iJ, and is therefore only valid for 
Z 3> 1 or at sufficiently late times. We show the resulting ec- 
centricity distribution for % = 0.5, eb = 0.3 and flb = 10 AU 
in Fig. [1] These results compare very well to the numerical 
results (see Fig. [5]). 

In the limit ef ^ \J Z — leg, some simple relations can 
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Figure 2. Equilibrium gas eccentricity distribution, analytically 
derived from Eq. for eb = 0.3, and qy, = 0.5 with /3 = -7/4 

(gray curve) or % = 0.234 with f3 = —1/2 (black curves). Also 
shown is the forced eccentricity Cf (black dotted line), which is 
the same for both cases. The dashed black curve shows the effect 
of moving the outer boundary inwards. All models use h = 0.05. 




be obtained: 



eo 



Z + 1 



tan{zuo) ~ 



1 ^ \/ "z^ tan(^A7g) 



(30) 
(31) 



where zug^ denotes the longitude of periastron of the gas 
disc. We see that in the limit of strong coupling to the 
gas {Z ^ 1), the planetesimals will tend to obtain the gas 
disc eccentricity. On the contrary, for very weak coupling 
(Z 1), their eccentricity will tend to zero. Note, however, 
that Eq. (O is not vaUd in the Umit Z ~* \\ from Eq. ((26)) 
it is in fact easy to see that eo et for Z ^ 1. 

An important conclusion is that, except in the special 
case Eg = et, friction with an eccentric gas disc always leads 
to some degree of differential orbital phasing between plan- 
etesimals of different sizes. This case is thus qualitatively not 
different from that of a circular gas disc. The magnitude of 
the differential orbital phasing effect will depend strongly on 
the gas eccentricity though. 



3.4 Gas disc eccentricity 

We now consider the crucial issue of the response of the gas 
disc itself to the perturbations. The equivalent of Eq. (|6]) for 



a gaseous disc reads: 
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where f3 = dlogE/dlogr, h — H/r, and primes denote dif- 
ferentiation with respect to r. It is possible to find equilib- 
rium solutions by setting the left-hand side to zero, after 
which the resulting second order ordinary differential equa- 
tion can be solved in terms of Bessel functions (for constant 
/3). Before doing so, we note that for the case /3 = —1/2, 
Eg — Bf is actually a solution (for appropriate boundary 
conditions) . 

It is not straightforward to supply realistic boundary 
conditions for this problem. It is conceivable that the gas 
disc eccentricity will tend to zero close to the star, and nu- 
merical hydrodynamical simulations (to which we intend to 
compare these results) will also force the eccentricity to be 
zero at the boundaries. Therefore we choose £'g(rin) — 0, 
where rin is the inner radius of the disc. The outer bound- 
ary poses a bigger problem, because the disc will be trun- 
cated. Since linear theory does not include the truncation 
of the disc, we vary the outer boundary condition between 
Eg{rout) = 0, where rout is the outer radius of the disc and 
£-g(''trunc) = 0, where rtrunc is the approximate truncation 
radius of the disc, obtained from hydrodynamical simula- 
tions. Note that since vj{rin) = ii7(rout) = 0, the disc will be 
aligned with the binary everywhere according to Eq. (|32|) . 
A possible twist in the disc does not change the results sig- 
nificantly. 

The results are shown in Fig. (2] for two binary config- 
urations: Qh — 0.5, eb — 0.3, P — —7/4 (gray curve), and 
qh = 0.234, Bb = 0.3, fH — —1/2 (black curves). Both cases 
will be studied numerically in Sect. [5] We put the bound- 
aries at rin ~ 0.025 ah and rout = 0.4 Ob, the same as will be 
used in Sect.[5]in the hydrodynamical simulations. For these 
binary parameters, the outer boundary is well beyond the 
truncation radius, while the inner boundary is well inside 
the planet-forming region (for ah ~ 10 AU). 

The steep density profile (denoted by the gray curve 
in Fig. (21) gives rise to an eccentricity well above ef for a 
large part of the disc. Reducing the inner radius of the disc 
by a factor of 2 has virtually no effect on the eccentricity. 
This is not true for the outer boundary, as we will see be- 
low. For P — —1/2, the eccentricity approaches ef far from 
the boundaries. This may be important, because in this case 
there would be no differential orbital phasing in the planet- 
forming region. However, this strongly depends on the lo- 
cation of the outer boundary. If we suppose that the disc 
is circular at the truncation radius, which, for this binary 
configuration, is located around rtrunc = 0.3 ah, the pic- 
ture changes drastically (see the dashed line in Fig. [2} . In 
this case, the disc remains nearly circular everywhere, and 
we will see in Sect.[5]that this closely resembles the solution 
obtained from numerical hydrodynam ical simulations. Also, 
there is a strong dependence on h, with larger eccentricities 
for thinner discs. 

However, there are three possibly important physical 
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Figure 3. Circular gas disc case, for a tight binary % = 0.5, 
e(j = 0.3 and a^, = 10 AU. Shown are analytical estimates of 
mean encounter velocities between planetesimals, at 1 AU from 
the primary, as a function of their size, derived using Eq. I I33I I. 
We consider target bodies of two different sizes, one with R = 
1 km and the other R = 5 km, and derive for each case collision 
velocities due to impacting obje cts in the 1-10 km size range. 
Overplotted are the results of ,Thebault et al-lliooel (their Table 
3). 

effects not included in Eq. (|32}. First of all, tidal effects are 
neglected, which means that the effect of disc truncation is 
not taken into account. Therefore /3 = /3(r) is unknown, and 
will change dramatically near the disc edge. 

There are also two possible eccentricity excitation mech- 
anisms not included in Eq. (I32p . There is the effect of the 
3:1 re sonance (|Lub ow 199 1,), and also t he viscous oversta- 
bility (|Katolll978l : iLatter fc Ogilviell2006l ') . Both mechanisms 
can add a significant free eccentricity to the disc. We rely 
on numerical hydrodynamical simulations to solve for the 
eccentricity of the gas disc. 

3.5 Encounter velocities 

A simple model, where one considers the collision between 
two kinds of bodies with the same semi-major axis in the 
guiding centre approximation, suggests that their mean en- 
counter velocity is proportional to — i72|- If we take 

JA^ = an^\Ei - E2\, (33) 

then this prediction reduces to the standard relation for bod- 
ies with mean eccentricity Cf and randomized periastra of 
iLissauer fc Stewart! l|l993h : 

M = ef. (34) 

When assuming a given eccentricity for the gas disc, 
Eq. ((33}, together with Eq. ((24}, can be used to predict 
the stationary solution for mutual encounter velocities. For 
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Figure 4. Planetesimal eccentricity (grey dots) and average en- 
counter velocity (solid line) distribution in a gas-free disc, for the 
tight binary case. 

the specific case of a circular gas disc, it gives a remark- 
ably good fit to the numerica l results obtain e d, on ce the 
stationary state is reached, bv lThebault et al] (|2006l l (Fig. 
(3}. Note however that this set of equations cannot predict 
if, and how fast the stationary state will be reached. It also 
cannot account for the possibly crucial effect of encounters 
due to high-eccentricity objects coming from other regions 
of the system (in case of orbital crossing). These equations 
are nevertheless very useful in providing a reference to which 
to compare numerical results obtained with an evolving gas 
disc. Such a comparison can provide information on whether 
the gas disc eccentricity dominates the coUisional evolution 
of the planetesimal population. 



4 TEST CASES: GAS FREE AND 
AXISYMMETRIC CASES 

We start by reproducing the results of lThebault et ahl (|2006l ) 
for the gas-free case and for an axisymmetric gas disc. These 
results can then serve as a reference for the full model, but 
they also provide good test cases for our method. 

4.1 Gas-free disc 

In the case of no gas drag, the eccentricity evolution of the 
planetesimals is governed by secular perturbations due to 
the binary only. These perturbations give rise to eccentric- 
ity oscillations, with a wavelength that decreases with time. 
Due to the strong orbital phasing large eccentricities may 
arise with relatively low encounter velocities. However, as 
time goes by the eccentricity oscillations become so narrow 

^ For the sake of comparison, we display here results for the 
same specific b i nary configuration considered as an example by 
iThebault et al.l ||2006|) . i.e., ab = 10 AU, % = 0.5 and eb = 0.3. 
This corresponds to a close and eccentric binary, for which the 
inner 1 — 2 AU region is highly perturbed by the companion star 



8 S.-J. 



Paardekooper, P. Thebault and G. Mellema 



a (AU) a (AU) 

0.6 0.8 1.0 1.2 1.4 1.6 1.8 0.6 0.8 1.0 1.2 1.4 1.6 1.8 




r (AU) t (yr) 

Figure 5. Planetesimal evolution in an axisymmetric gas disc, with for the same binary parameters as in Fig. |2] (see also the 

top left panel). Top left: longitude of periastron distribution after 5680 yr. Top right: eccentricity distribution after 5680 yr, with the black 
line indicating the forced eccentricity ef. Bottom left: distribution of encounter velocities for 3 different impactor-target planetesimal 
pairs: 2 corresponding to equal sized impacting objects and the third one to a 1 km impactor hitting a 5 km target. Also shown is the 
limiting maximum {dv) for accreting encounters between 1 and 5 km bodies (dotted line). Bottom right: encounter velocities at 1 AU as 
a function of time for the same impactor-target pairs. 



that orbital crossing between high and low eccetricity bod- 
ies eventually occur, resulting in a sudden and dramatic in- 
crease of encounter velocities. Orbital crossing first occurs 
in the outer regions (closer to the companion star), the ra- 
dial location at which orbital crossing occurs progressively 
moving inward with time. The time at which orbital cross- 
ing occurs at a given location can be calculated analyticall y 
following the empirical derivation of iThebault et ah! (|2006l ). 

In Fig. |4]we show the eccentricity distribution of plan- 
etesimals after 6200 yr (200 binary orbits) in a tight binary 
system with parameters as indicated above the figure. The 
eccentricity oscillation can be identified easily. In the outer 
part of the disc, several particles have been captured into res- 



onances. Note that we only consider here the inner ^ 2 AU 
part of the disc, since most planetesimal orbits beyond this 
distance are dynamically unstable. 



The black solid line in Fig. |4] indicates the mean en- 
counter velocities as a function of distance. The radius where 
orbital crossing occurs can be found near 1.1 AU (note the 
very sharp transition from (Av) ~ to very high veloc- 
ity values in less than 0.1 AU around the crossing location. 
These results quantit atively agree with the calculations of 
iThebault et"alT (|2006l ) ( see their Fig. 2). 
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Figure 6. Surface density of the gas disc, with Eg oc r~^/* after 50 binary orbits (fs 1600 yr), for the same binary parameters as in 
Fig- HI Left panel: two-dimensional distribution of logjQ S. The binary companion is at apo-astron, at (x,y) = (—13,0). Right panel: 
azimuthally averaged surface density. The density scale is arbitrary, for comparison the initial condition is shown. 



4.2 Axisymmetric gas disc 

iThebault et al.l ()2006l ) included effects of gas drag from an 
axisymmetric gas disc. In this secti on, we use the same gas 
disc model as IThebault et ahl l|2006l l and do not let it evolve. 
We divide the 10000 planetesimals into two groups of 1 km 
and 5 km in size. We consider the same tight binary con- 
figuration as in the gas free case: Ob = 10 , eb = 0.3 and 
Qh = 0.5. In Fig. O we show the resulting orbital and col- 
lisional evolution, in the ^ 2 AU region corresponding to 
stable orbits. These c an be compared to Figs. 6 and 7 of 
IThebault et"ai] (|2006l l. 

In the top panels of Fi g. [5l we observe the w ell known 
fe.g. lMarzari fc Sch oh 2000.: IThebault et al.ll2004l ) double ef- 
fect of gas drag on planetesimal orbits, i.e., damping of the 
eccentricities and periastron alignment. At first sight, this 
may seem advantageous for planetesimal accretion. How- 
ever, the equilibrium between eccentricity excitation and 
damping depends on the size of the planetesimals. Not only 
are the eccentricities of the smaller bodies damped more 
readily than for the larger bodies, they are damped towards 
a different equilibrium. This is also true for the equilibrium 
value of the periastron alignment. This differential phasing 
is responsible fo r high collision veloci ties between particles 
of different sizes (|Thebault et al.ll2006l l. Note also that in the 
inner disc, where the eccentricity oscillations are damped by 
gas friction, the eccentricities and longitudes of periastron 
of the planetesimals agree with the analytical results in Fig. 

m. 

These consequences of gaseous friction are clearly illus- 
trated in the bottom panels of Fig. [5] which shows the evo- 



lution of the average encounter velocities (Av) ri^r2, for the 
cases -Ri = -R2 = 1 km, Ri = R2 = 5 km, _Ri = 1 km and 
R2 = 5 km. A first obvious effect, clearly seen for the case 
of equal-sized objects, is that gas drag works against orbital 
crossing: stronger gas friction moves the radius at which or- 
bital crossing occurs outward. At 1 AU, no orbital crossing 
occurs after 3000 yr for both planetesimal sizes and rela- 
tive velocities are still very low (see the bottom left panel 
of Fig. [5)1. In sharp contrast, encounter velocities between 
bodies of different size are high throughout the disc. In the 
bottom right panel of Fig. [S] we see that the encounter ve- 
locities reach a stationary value after approximately 5000 
years. T he equilibrium is in pe rfect agreement with the re- 
sult from IThebault et al.l (|2006l l. In the bottom left panel of 
Fig. [5] we also show the limiting value of (dv) for accreting 
collisions between planetesimals of 1 and 5 km (dotted line). 
This is a conservative estimate, b ased on the most optimistic 
value for accretion considered bv IThebault" et al. I l|2006l l. 



5 RESULTS, FULL MODEL 
5.1 Gas disc evolution 

We start the discussion of the full model by considering 
the gas evolution alone. In Fig. |B]we show the surface den- 
sity of the gas disc after 50 binary orbits using the min- 
mod flux limiter. The disc is truncated by the binary at ap- 
pr oximately the radiu s pred icted by the empirical formula 
of iHolman fc WiegertI (ligogt ). A large part (approximately 
20%) of the material is expelled from the system, but from 
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Figure 7. Gas disc eccentricity (top panel) and longitude of pe- 
riastron (bottom panel) for the excited state, obtained with the 
soft superbee flux limiter (solid line) and the quiet state, ob- 
tained with the minmod flux limiter (dotted line). The binary 
parameters correspond here to the wider binary case of 7 Cephei: 
gb = 0.234, ab = 20 AU, ej, = 0.3, and the initial surface density 
has a shallow radial dependence Eo DC r-1/2. 

The thick solid line 
indicates the forced eccentricity derived from Eq. II12I I. 



the right panel of Fig. |6]it is also clear that the disc is com- 
pressed: the density at 1 AU is increased by approximately 
a factor of 3. This means that in order to compare with the 
results of the previous section, we should rescale the density 
in the full model to obtain the same magnitude of the drag 
force at 1 AU. 

The second thing that is apparent from Fig.|5]is the ap- 
pearance of spiral density waves. The morph ology is strongly 
depe ndent on the phase of the binary (see iKlev fc NelsorJ 
I2OO7I ). with strong shocks appearing just after periastron. 
The effect of these waves in a circular binary on the or- 
bital elements of planetesimals was recently studied in 
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Figure 8. Evolution of the mean disc eccentricity and longitude 
of periatron, for the wide binary case q^, = 0.234, = 0.3 and 
= 20 AU, with Sq oc r~^^^ , using two different flux limiters. 



ICiecielag et"al] l|2007h . where it was shown that only the 
smallest bodies {R « 100 m) are affected. 

Finally, it is clear from the left panel of Fig. [6] that the 
gas disc becomes eccentric. This is always true, but tests 
show that the two considered flux limiters give very differ- 
ent results. The amplitude of the eccentricity increase de- 
pends dramatically on the amount of wave damping, as is 
illustrated in Fig. [T] where we compare the eccentricity dis- 
tribution after 50 binary orbits for the 7 Cephei binary con- 
figuration, (5b = 0.234, Ob — 20 AU, eb = 0.3) for results 
obtained with our two different flux limiters. 

The diffusive minmod limiter gives rise to a low eccen- 
tricity (gg ^ 0.05 and even ^ 0.02 in r < 4 AU region), 
stationary disc state that we henceforth call the quiet state. 
It is a robust state, numerically, to changes in resolution and 
boundary conditions (reflecting, non-reflecting). In fact, this 
solution approximately obeys the time-independent version 
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of Eq. (|32|l . as we will show below. The maximum eccen- 
tricity does depend on the location of the inner boundary, 
however, since the boundary conditions enforce the eccen- 
tricity to be zero at the boundary. We have found that the 
location of the inner boundary should be at least as far in 
as 0.025 flb to get converged results. 

The soft superbee flux limiter gives rise to a different 
disc state, with a large free eccentricity, reaching values up 
to eg ~ 0.2. We will refer to this state as th e excited state. 
This state is similar to the one discussed in iKlev fc NelsonI 
l|2007h . Its exact behavior depends strongly on resolution 
and boundary conditions, but the overall picture is the same: 
high eccentricity, and the disc starts to precess. This is il- 
lustrated in Fig. [S] where we show the mass-averaged eccen- 
tricity and longitude of periastron for the disc. Within 20 
orbits, the disc obtains a large eccentricity, after which the 
disc starts to precess slowly in a retrograde fashion. The ec- 
centricity oscillates with the same period as the precession. 
The retrograde precession indica tes that pressure forces are 
dominant (see iPataaloizo ul l2005h . In Fig. [H the quiet disc 
quickly reaches a steady eccentricity distribution. Note that 
due to the twist in the disc (see Fig.[7|, which is dynamically 
unimportant, the average longitude of periastron approaches 
neither nor tt, but somewhere in between. 

Note that for both states the disc eccentricity does not 
equal the forced eccentricity, denoted by the thick solid line 
in Fig. [7] For the quiet state, we have checked that the gas 
eccentricity approximately follows Eq. (|32|l . The fact that, 
despite the density being proportional to r'^''^, the gas ec- 
centricity does not tend to the forced eccentricity is due to 
the behavior near the truncation radius. This is basically 
where the boundary condition for Eq. (|32|) is determined, 
which, as shown in Fig. [2] can easily lead to a strongly re- 
duced disc eccentricity. According to the analytical study 
of Sect. |3l this means that we can expect large encounter 
velocities between planetesimals in both cases. 

The reason for the large differences between the runs 
with different flux limiters is that the mechanism responsible 
for generating a large free eccentricity involves an eccentric 
resonance near the out er edg e of the disc. Basically, it is the 
mechanism outli ned in iLub ow (1991), but for an eccentric 
companion orbit. iHeemskerkl (.1994i ) showed that eccentricity 
excitation through this instability strongly depends on the 
way the disc is truncated, which in turn depends strongly on 
wave damping. The two different limiters give different wave 
damping, hence our results, see Appendix A. Only for small 
wave damping does the disc obtain a large free eccentricity. 
It is expected that this also depends on the magnitude of 
the physical viscosity, since this affects the shape of the disc 
edge. We provide some additional comments in Sect.[SJ but 
leave a detailed discussion of this problem to a future pub- 
lication. We now proceed to analyze the influence of both 
disc states on planetesimal accretion. 

5.2 Planetesimal evolution 

5.2.1 Quiet disc case 

We first consider the quiet disc case, for which the situation 
should a priori be the closest to the non-evolving axisym- 
metric case. 

Fig. [9] shows the evolution of planetesimal orbital pa- 



rameters and encounter velocities for the same tight bi- 
nary parameters as in Fig. [5] Comparing the top panels of 
Figs. [S] and [5] we see that, as expected, the largest differ- 
ences occur in the inner regions of the disc, typically within 
~ 0.8 AU, where gas drag effects are the most important. 
The 1-km planetesimals approach the gas eccentricity to- 
wards the inner boundary, as expected. Interestingly, the 
5-km bodies show a large jump in longitude of periastron 
around r = 0.7 AU, accompanied by a drop in eccentricity. 
From the top left panel of Fig. [9] we see that this happens 
where the longitude of periastron of the gas disc amounts 
to 7r/2. Around this location, depending on the value of Z, 
the denominator of Eq. (|27|) will approach zero, causing a 
large jump in zu. For the 5 km planetesimals, Z ~ 3 at this 
location, which causes a drop in eccentricity (see Eq. (|26p ). 
This is not true for the 1 km planetesimals, and therefore 
there is a large eccentricity difference around r = 0.7 AU. 

In terms of encounter velocities, these different be- 
haviours of 1 and 5 km bodies in the innermost regions logi- 
cally translate into higher Aw than in the axisymmetric case 
(see the bottom-left panel of Figs. [5] and [9]). At 1 AU, the 
equilibrium encounter velocities are approximately a factor 
of 2 higher than for the case of a circular gas disc. However, 
in the outer disc, beyond ~1 AU, differences with the ax- 
isymmetric case are much smaller. Although the gas eccen- 
tricity is higher in these regions, the gas density is not high 
enough to significantly affect the behaviour of km-sized plan- 
etesimals. Beyond r = 1.4 AU, the dynamical evolution of 
the planetesimal population becomes indistinguishable from 
the circular gas disc case. 

We now turn our attention to binary parameters that 
match those of 7 Ceph: gb = 0.234, ab = 20 AU and = 0.3. 
From Fig. [7] we see that the gas disc eccentricity is very 
small, eg < 0.02 almost everywhere. In the top panels of 
Fig. [To] we show the distribution of longitude of periastron 
and eccentricity after 10* yr, when the system has reached a 
steady state. We find that in the whole r ^ 0.8AU region, the 
equilibrium (Au) for collisions between 1 and 5 km objects is 
^ 200 ms~^. This is still high enough to correspond to erod- 
ing impacts for all tes ted collision outcome prescriptions of 
iThebault e'ta\\ l|200d l (see bottom-left panel of FiglTUI. Di- 
rect comparison with Fig. [5] is here difficult, because the 
binary parameters are different. We thus performed an ad- 
ditional axisymmetric gas disc test simulation which showed 
that the encounter velocities are the same, within 10%, as 
for the present quiet state run. This is an indication that 
the spiral waves, that do extend all the way in, are indeed of 
minor importance regarding impact velocities. The circular 
gas disc case is then a relatively good approximation. 

5.2.2 Excited disc case 

In Fig. [Tl]we show the results for the excited disc state, for 
the 7 Cephei like binary. It is immediately clear that, even 
for this wide binary case, the planetesimals react rather vi- 
olently to the large eccentricity of the gas disc. Inside 1 AU, 
the 1 km planetesimals are dragged along with the gas and 
end up on highly eccentric orbits matching that of the gas. 
The 5 km planetesimals follow the same trend, but are more 
loosely coupled to the gas streamlines and their eccentrici- 
ties never match that of the gas. For the periastra, however, 
we observe a quasi-perfect alignment with that of the gas for 
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Figure 9. Planetesimal evolution in an evolving gas disc (the quiet state), with a steep radial gas profile in Eg oc r~'^ 1^, for the same 
tight binary parameters as in Fig. |5] Top left; longitude of periastron distribution after 5680 yr. Top right: eccentricity distribution after 
5680 yr. Bottom left; distribution of encounter velocities for 3 different impactor-target planetesimal pairs; 2 corresponding to equal sized 
impacting objects and the third one to a 1 km impactor hitting a 5 km target. Also shown is the limiting {dv) for accreting encounters 
(dotted line). Bottom right; encounter velocities at 1 AU as a function of time for the same 3 impactor-target pairs. 



both planetesimal sizes. Note that this rather abruptly hap- 
pens where V3^ « 37r/2, which makes | tan(tj7g)| very large. 
From Eq. l(3T]l we see that in this case, for which we indeed 
have eg ^ Cf, we expect all planetesimals to align with the 
gas, independently of Z. Outside 1 AU, 1 km planetesimals 
start to decouple from the gas, while the bigger bodies do so 
already at ~ 0.7 AU. In the whole 1 ^ r ^ 3 AU region, dif- 
ferences between the equilibrium eccentricities for both sizes 
are very large, culminating at ~ 2 AU, where eikm ~ 2 eskm- 

The differences in E between the two particle sizes in 
the planet-forming region around 1 AU lead to very high 
encounter velocities, exceeding 500 ms^^ almost everywhere 
(bottom panel of Fig. fTTj) . Note that due to the precession of 



the disc, is time dependent, as appears clearly from the 
oscillations in Fig. [8l However, this precession timescale is 
shorter than the time it takes for the planetesimals to settle 
into their equilibrium distribution. Therefore, the planetes- 
imals will feel an average gas eccentricity and eventually 
settle into a steady eccentricity and (Au) distribution. Also, 
the disc at 1 AU does not really participate in the precession, 
possibly due to the fact that it is close to the inner bound- 
ary (see also below) . This is also apparent from Fig. [7] For 
comparison, we also show the encounter velocities evolution 
at 2 AU in the bottom right panel of Fig. lTTI We see that, at 
this larger distance, the encounter velocities oscillate, with 
a period which is that of the gas disc precession (see Fig.jS]), 
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Figure 10. Planetesimal evolution in an evolving gas disc (the quiet state), with a shallower radial gas profile Sq oc r~^/^, and for 
the wider, 7 Cephei like binary case. Top left: longitude of periastron distribution. Top right: eccentricity distribution. Bottom left: 
distribution of encounter velocities for 3 different impactor-target planetesimal pairs: 2 corresponding to equal sized impacting objects 
and the third one to a 1 km impactor hitting a 5 km target. Also shown is the limiting (dv) for accreting encounters (dotted line). 
Bottom right: encounter velocities at 1 AU as a function of time for the same 3 impactor-target pairs. In addition, we show the encounter 
velocities at 2 AU for a 1 km impactor hitting a 5 km target (dotted line). 



but after approximately lO^yrs a meaningful average can be 
defined. We note that, at a given location in the disc, en- 
counter velocities are close to those predicted by Eq. (|33|l . 
when plugging in the average value of Eg at this location (see 
also Fig. I12|l . This again shows that it is the eccentricity of 
the gas disc that governs the impact velocities. 

The bottom left panel of Fig. [11] shows that encounter 
velocities in the inner regions, which are due to (gas drag 
induced) differential orbital phasing, come close to those in 
the outer disc, which are due to pure gravitational orbital 
crossing. For lkm-5km pairs, impact velocities are on aver- 
age 5 times higher than those in the circular disc run. They 



are thus far above the upper limit for accreting impacts for 
the whole simulated region 0.6 < r < 6 AU. 

We have also run the soft superbee flux limiter for 
the tight binary case. Differences with the circular case are 
slightly larger than for the quiet state, but the gas disc is not 
able to obtain a large free eccentricity in this case, for two 
reasons. First, numerical and analytical experiments suggest 
that the maximum eccentricity (free or forced) that the gas 
disc can reach is largely determined by the boundaries (see 
also Fig. m . The quiet state of the tight binary comes close 
to obtaining this maximum of eg ~ 0.3, and therefore the 
excited state is not very much different in this case. Sec- 
ond, because the gas disc is truncated closer to the primary, 
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Figure 11. Planetesimal evolution in an evolving gas disc (the excited state), with for the same binary parameters as in 

Fig. 1101 Top left: longitude of periastron distribution. Top right: eccentricity distribution. Bottom left: distribution of encounter velocities 
for 3 different impactor-target planetesimal pairs: 2 corresponding to equal sized impacting objects and the third one to a 1 km impactor 
hitting a 5 km target. Also shown is the limiting (dv) for accreting encounters (dotted line). Bottom right: encounter velocities at 1 AU 
as a function of time for the same 3 impactor-target pairs. In addition, we show the encounter velocities at 2 AU for a 1 km impactor 
hitting a 5 km target (dotted line). 



due to the stronger perturbations from the tight binary, the 
mechanism for generating a large free eccentricity is reduced 
in strength compared to the wide binary case. Therefore, en- 
counter velocities do not vary more than 50 % between the 
quiet and excited state for the tight binary, meaning that 
velocities are here a factor of ~ 3 higher than in the circular 
gas disc case. 

As already mentioned, numerical experiments have 
shown that the position of the inner boundary determines 
up to what radius the gas disc obtains its large eccentricity. 
For example, in Fig. [71 the eccentricity drop inward of 2 AU 
is directly connected to the location of the inner boundary 
at 0.25 AU. Exploring this parameter is outside the scope 



of the present paper, since moving the inner boundary fur- 
ther inwards is computationally very expensive. However, it 
is important to note that it would only act to further in- 
crease the encounter velocities. More work is necessary to 
study the behavior of the excited disc state close to the in- 
ner edge of the disc, numerical or physical. Here, we focus 
on the influence of a gas disc with a large free eccentricity 
on planetesimal collisions only, and we have shown that it 
is disastrous for planetesimal growth. 
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Figure 12. Encounter velocities, calculated using Eq. I I33I I for 
the mean eccentricity of the excited gas disc state, for binary 
parameters = 0.234, = 20 AU and = 0.3. Black lines 
were calculated using _Eg = 0.05, appropriate for 1 AU (see Fig. 
0, grey lines using i?g = 0.15, appropriate for 2 AU. 



6 DISCUSSION 

6.1 The impact of gas disc eccentricity on 
encounter velocities 

Our results indicate that planetesimal encounter velocities, 
which result from differential orbital phasing due to gas 
drag, critically depend on the eccentricity of the gas disc. Al- 
though there is a theoretical possibility for differential phas- 
ing to be avoided, when gas streamlines follow the forced 
dynamical orbits, we have shown that this situation does not 
occur in practice. First of all, Eq. (|32|l shows that we can 
only have — et when the gas density follows E oc r^^^^ . 
Even if this happens to be the case, there is a strong depen- 
dence on what happens at the inner and outer boundaries of 
the disc (see Fig. [2} , in combination with the disc thickness. 
More importantly, this = a condition was never encoun- 
tered in our numerical hydrodynamical simulations. Since 
we could not explore all free parameters, we cannot rule 
out the possibility of obtaining _Bg = ef very locally when 
tuning in the right disc parameters, but this would have to 
be a fortuitous coincidence. Even more, if the disc reaches 
a large free eccentricity, then it will necessarily precess, so 
that iJg = ef would not be a steady state. 

Not only does this Eg — ef state never occur, but for all 
4 cases here explored (2 flux limiters -I- 2 binaries configura- 
tions), we always obtain, at all radial distances, (Aii) equi- 
librium values higher than for a static circular gas disc case. 
For the quiet state, differences with the static case remain 
relatively limited: velocity distributions are almost identi- 
cal for a binary with 20 AU separation and are enhanced 
by slightly less than a factor of 2 for a tight ah — 10 AU 
companion. For the excited state, however, differences with 
the static disc case can become large. They already exceed 
a factor of 2 for the wide (20 AU) binary case and reach a 



factor of 3 for the tight binary. Interestingly, even for the 
case where eg turns out be relatively close to et, i.e., in the 
inner disc in the quiet state/tight binary run (see Fig. O, 
we still observe an important differential phasing between 
particles of different sizes and thus high (Av) values. This is 
entirely due to the fact that the gas disc is not aligned with 
the binary in that region. This is a further indication that 
the parameter phase space around the exact iJg = et condi- 
tion where differential phasing is weaker than for a circular 
disc is very limited. 

The computational cost of the presented numerical sim- 
ulatio ns inhibits an extended parameter study similar to the 
one m iThebault et"al] (|2006h for the static circular gas disc 
case. As already mentioned, we had to restrict ourselves to 
two representative binary configurations as well as to two 
planetesimal sizes, 1 and 5 km. However, we have checked 
(see Fig |12|l that for all numerically explored cases, Eq. (|33p 
gives good estimates of equilibrium encounter velocities once 
given the value for the mean gas disc eccentricity iJg (which 
is obtained from Fig. [7] and an equivalent plot for the quiet 
state). As a first approximation, this equation can thus be 
used as a tool to estimate (An) between planetesimal of any 
sizes -Ri and R2, for any given binary orbital configuration. 
The only unknown input parameter is Eg at a given location 
in the disc, but this value can be retrieved from independent 
pure hydro-simulations, which are less time consuming than 
full gas-fplanetesimal runs. 

For the present binary configurations, we use Eq. (|33|l 
to estimate impact velocities for planetesimals bigger than 
the ones explored in the runs (see Fig |12l for the excited 
case). This could be of interest for two reasons: 

(i) Sizes in the "initial" planetesimal population, however 
fuzzy this concept might be, are very poorly constrained by 
planet-formation scenarios. 

(ii) The steady state with equilibrium {Av)ri^r2 is not 
reached instantly and might take a few 10^ yrs to settle (see 
Figs. O and lll|l , which might leave time for some planetesi- 
mals to start to grow and reach bigger sizes. 

(for more detailed discussions on these 2 issues, see 
IThebault eral] |2006| ). From Fig. [12l we see that, for the 
7 Cephei-like binary, impact velocities quickly rise to very 
high values for larger bodies. Even though exact collision 
outcomes for such large objects remain poorly constrained, 
such extremely high velocities are very likely to lead to 
mass erosion or shattering of the impactors. As an example, 
for the same binary configuration (except a slighly higher 
Qh value) and for col lisions between 20 and 50 km bodies, 
IThebault et"al] {200& \ find {Av) ~ 300 m s with a circular 
gas disc, a value for which no clear conclusions in terms of 
accreting vs. erosive impacts can be drawn (see Fig. 9 of that 
paper). For the present excited case, however, values larger 
than 10"^ ms^^ are reached, which are without doubt high 
enough to lead to net mass loss after the impact. 

6.2 The origin of gas disc eccentricity 

Additional test simulations with a circular binary orbit 
showed similar behavior as for the eccentric binary case. 
Depending on the amount of wave damping, the disc can 
quickly become eccentric. We have check ed that in thi s case 
it is the 3:1 eccentric Lindblad resonance (|Lubowlll99il ) that 
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drives the eccentricity. Removing the m = 3 component from 
the companion potential resulted in almost circular discs for 
all flux limiters. This indicates that a similar mechanism 
operates in the case of an eccentric binary. When the disc 
does become eccentric, the circular binary case behaves in a 
similar fashion as the eccentric binary case. This is because 
l^'gl S> efi the eccentricity evolution of the planetesimals is 
dominated by the free eccentricity of the gas disc. Whether 
the 3:1 resonance can induce a large eccentricity in the gas 
disc strongly depends on the truncation of th e disc. It is no t 
expected to happen for equal- mass binaries (|Lubowlll99l} ). 
which was the case studied bv ICiecielag et aL I (|2007l)" for a 
non-evolving gas disc. More work is necessary to find out at 
which mass ratio this instability sets in and what values of 
eg can be reached. 

We have not found evidence for a viscous overstabil- 
ity operating in the disc. However, there exists a complex 
interplay between viscosity, tidal effects and eccentricity ex- 
citation. If the disc is truncated inside the main eccentricity- 
generating resonance, the disc will not obtain a large free ec- 
centricity. However, a strong viscosity may spread the outer 
edge of the disc into the resonance, and this complex prob- 
lem clearly deserves more study. Our formalism, in particu- 
lar Eq. (|33p . allows for a direct link between gas disc eccen- 
tricity and planetesimal encounter velocities, at least in the 
region where orbital crossing can be neglected. This may 
ease the computational cost of future studies, as simula- 
tions with gas only are required to investigate planetesimal 
encounters. 

An important parameter governing gas friction is the 
ambient gas density. In our simulations, we always rescale 
the density at 1 AU to pg — 1. 4 ■ 10~^ g cm~ ' ^, in order 
to compare with previous work (|Thebault et al.ll2006l ). the 
actual density may be very different. Even for single stars 
the disc mass is relatively poorly constrained, and the effect 
of the binary further complicates things. As the gas density 
only appears in the drag force parameter K, which also con- 
tains the planetesimal size, our results can always be scaled 
to different gas densities by considering different planetesi- 
mal sizes. 

The effects of self-gravity are usually neglected in cal- 
culations of protoplanetary discs, because unless the disc 
is very massive (Md ~ 0.1 Mq) its dynamical influence is 
small. The parameter measuring t he imp ortance of self- 
gravity is the Toomre Q parameter (|Toomrail964i l : 



= 



ttGE' 



(35) 



where k is the epicyclic frequency. For a locally isothermal 
Keplerian disc with a constant aspect ratio we have: 



Q = 



h_ 

where 



Mr: 



(36) 



(37) 



is a measure of the disc mass. Taking a disc mass of the order 
of one Jupiter mass we find Q ~ 50, which indicates that 
self-gravity is not important. However, Toomre's stability 
criterion is based on short wavelengths. For global modes, 
self-gravity can d ominate over pressure much more easily 
l|Papaloizoull2002h . Basically, the parameter measuring the 



importance of self-gravity becomes /iQ, which is of order 
unity for our model disc. Therefore, future models should 
include self-gravity of the gas. 

Within our analytical framework, taking into account 
gravitational forces due to the gas amounts to adapting the 
forcing potentials appearing in Eq. (|6]) to include perturba- 
tions due to the gas. These will involve integrals of the gas 
surface density over the disc, where for $2 the integrand 
will depend on Eg. For large enough disc mass and disc 
eccentricity, the forcing due to the eccentric gas disc will 
dominate over the binary forcing. We may expect this to be 
the case in the excited state. Although a detailed analysis 
of this case is beyond the scope of this paper, we comment 
that we do not expect this case to differ qualitatively from 
the case studied here. Again, there will be an equilibrium 
eccentricity distribution for the gas, which will be harder to 
compute because now an integro-differential equation has 
to be solved, and an equilibrium eccentricity distribution 
for the planetesimals. In general, these will be different due 
to the spatial derivatives appearing in the gas eccentricity 
equation, which will give rise to differential orbital phasing. 
Again, there may be a special case for which Eg = Eq for 
all planetesimal sizes, but it is not clear whether this state 
is reached more easily than for the cases studied here. This 
clearly deserves more study. 



7 SUMMARY AND CONCLUSION 

We present the first simulations of planetesimal dynamics in 
binary systems including full gas disc dynamics. As in past 
studies with static axisymmetric gas discs, we confirm the 
crucial role of differential orbital phasing due to gas drag in 
inducing large encounter velocities between bodies of differ- 
ent sizes. Interestingly, while there is a theoretical possibil- 
ity for differential phasing to be less pronounced than in the 
axisymmetric case (if -Ej ~ e/), our numerical exploration 
shows that this case never occurs in practice: we always find 
stronger differential phasing than for a simplified static cir- 
cular gas disc. 

The level of the phasing is connected to the eccentric- 
ity reached by the gas disc. Depending on the amount of 
wave damping, the disc either enters a nearly steady state, 
for which encounter velocities are within a factor of 2 differ- 
ence from the case of a circular gas disc, or the disc enters a 
highly eccentric state with strong precession, for which the 
encounter velocities can go up by almost an order of magni- 
tude in the most extreme cases. It is important to point out 
that, for both cases, the global qualitative effect is always 
to increase {Aw){fli_fl2) for _Ri 7^ R2. 

Taking into account the gas disc's evolution thus leads 
to a dynamical environment which is less favourable to plan- 
etesimal accretion. In this respect, despite of the fact that 
the high CPU cost of these simulations prevented us from 
performing a thourough exploration of all binary parame- 
ters, it is likely that the (ab,eb) space for which planetesimal 
accretion is possible is more reduced than w as found in pre- 
vious studies with an axisymmetric disc (e.g. lThebault et al] 
I2OO6I) . 
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APPENDIX A: FLUX LIMITERS AND 
ECCENTRIC DISCS 

In view of the large influence of the choice of flux limiter on 
the resulting disc structure in a binary system, it is appro- 
priate to review their use in numer i cal hy drodynamics. Our 
discussion closely follows iLeveaud (|l990l ). For most prob- 
lems the choice of flux limiter does not play a major role 
and thus it is not common practice to check the effect of the 
choice of flux limiter. Our results show that it would be wise 
to always do so, especially in highly non-linear problems. 

Al Linear advection equation 

We keep the discussion simple by considering the linear ad- 
vection equation, in stead of the non-linear system of Euler 
equations: 

du du ^ , . , , 

where u is the conserved quantity (eg. mass, momentum, 
energy) and a is the advection velocity, which we take to be 
a positive constant. We will denote the numerical solution, 
discretized on an equidistant grid {xi}, at time step n as 
{(7"}. An important difBculty arises in numerically solving 
this type of equation, because it allows for weak, or discon- 
tiimous, solutions. The infinite gradients at discontinuities 
are problematic for numerical difference schemes, leading to 
strong oscilla tions in the num erical solution. In fact, a the- 
orem due to iGodunovl (|l959l ) states that a linear method 
that preserves monoticity (i.e. does not lead to oscillations) 
can be at most first-order accurate. Therefore, while in re- 
gions of smooth fiow it is desirable to have a scheme that is 
second-order, in the presence of discontinuities or shocks it 
is necessary to switch to a first-order scheme. This is where 
the flux limiter enters the stage. 

Suppose we have a method for determining the numeri- 
cal flux F across a cell boundary. Moreover, we should have 
a second-order flux, F2 , and a flrst-order flux, Fi . In regions 
of smooth flow, we want to use F2 for higher accuracy, while 
near a discontinuity we want to use Fi to avoid spurious 
oscillations. To achieve this, we write the total flux as: 



F = Fi+${9) (F2-F1), 



(A2) 



where ^{0) is the flux limiter (not to be confused with the 
potential $ in the main text). In regions of smooth flow, 
"1> should be close to 1, while near discontinuities we want 
<I> = to switch to the flrst-order solution. Given a measure 
of the smoothness 0{U), it is possible to obtain bounds on 
the function al form of ^( 9) to prevent spurious oscillations 
near shocks (|Swebvlll98j ). Within these bounds, one is free 
to choose the flux limiter, and for all choices one obtains a 
method that is second-order in regions of smooth flow, and 
does not introduce oscillations near shocks. For Eq. (|A1|I . a 
natural choice for 9 is the ratio of two consecutive gradients: 



(A3) 



Note, however, that this measure of smoothness breaks down 
near extreme points of U. 

A different, and more graphical, interpretation of <1> is 
to regard it as a slope limiter. To evolve U for a time At, 
we can regard [/ as a piecewise constant function, shift it 
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by an amount of aAt, and finally integrate to obtain new 
cell averages. This leads to the upwind method (note that 
we assume a > 0): 



1.2 



ur 



(A4) 



where e = aAt /Ax is the Courant number. For numerical 
stability, we should have e < 1. The upwind scheme is a 
first-order method, and does not produce oscillations near 
discontinuities. 

We can improve on the upwind scheme by releasing the 
assumption of U to be piecewise constant, and in stead allow 
U to vary linearly within one grid cell with slope a. Applying 
the same technique (shifting and integrating this new U) 
leads to a different update for U: 

i)-^^2;e(l-e) (a^ - a"_i) . 



(A5) 



For a > 0, if we make the natural choice at — (f/i+i — 
Ui)/Ax, Eq. HA5|) reduces to the Lax-Wendroff method, 
which is second-order but leads to oscillations near shocks. 
If we in stead use 

+ 1 - Ur 



(A6) 



we end up with a numerical fiux of the form of Eq. (|A2[l . in 
which $ can now be interpreted as a slope limiter. In Fig. lAll 
we show the difference between three different slope limiters 
near a sharp gradient. The Lax-Wendroff slope (dotted lines) 
is obtained by setting $ = 1. The minmod and superbee 
limiters are examples of a general class of limiters of the 
form 



$(61) = max(0, min(l,p6'), min(0,p)). 



(A7) 



where for p = 1 the limiter is called minmod, and for p = 2 
we have the superbee limiter. In the main text, we have used 
p = 1.5, which we call the soft superbee limiter. For 1 ^ p ^ 
2, it can be shown that this limiter prevents oscillations in 
the numerical solution near discontinuities. Away from the 
sharp gradient, all limiters produce the same result. 

From Fig. lAll it is clear that the Lax-Wendroff slope 
will introduce oscillations near x = 0.45. The other two 
limiters will keep the solution monotonic, but with differ- 
ent slopes. In general, the superbee limiter will give rise to 
the steepest slopes possible without introducing oscillations 
in the solution. This way, numerical smearing of shocks is 
minimal, but at the same time, shallow gradients may be 
steepened artificially. In practice, it is safer to use a "softer" 
limiter, where one would expect that minmod will give re- 
sults that come closest to those obtained with non-shock- 
capturing schemes, which need explicit artificial viscosity to 
stabilize shocks. 



A2 Application to eccentric discs 

We expect the infiuence of the fiux limiter to be largest 
near discontinuities, since all limiters give $ = 1 in re- 
gions of smooth fiow. The problem at hand introduces strong 
non-linearity in the disc response. The spiral waves induced 
by the binary companion are in fact spiral shocks, and 
their strength is d irectly related to the growth of eccen- 
tricity in the disc (|Lubowl [l99ll ). Moreover, their dissipa- 
tion is related to the loc ation at which the disc is truncated 
l|Lin fc Papaloizoulll986l ), which also affects the strength of 
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Figure Al. Numerical solution U of Eq. IIAll l at time t„, and 
Ax = 0.05 (dots), together with the resulting interpolation for 
three different slope limiters. 

the instability (|Heemskerkll 199j ) . Finally, when eg > h, su- 
personic radial velocities occur. Since usually eg ^ near 
the boundaries of the computational domain, this will intro- 
duce shocks in the problem as well. 

It is clear, therefore, that shock dissipation plays a cru- 
cial role in the growth as well as in the damping of eccen- 
tricity in the gas disc. A diffusive flux limiter will reduce the 
effectiveness of eccentricity excitation, and will damp the 
obtained eccentricity more readily. Our results in the main 
text indicate that this leads to the mechanism involving the 
3:1 Lindblad resonance being unable to induce eccentricity 
in the disc when the minmod limiter is used. This strong in- 
fluence of the flux limiter is remarkable, and indicates that 
strongly non-linear effects play a large role in this prob- 
lem. Also, it is interesting that non-shock-capturing meth- 
ods seem to reproduce the results obtained with the soft 
superbee flux li miter rather than tho se obtained with the 
minmod limiter jKlev fc Nelsonll2007l ). The limiter used by 
iKlev fc Nelsoiil (|2007l ) is the so-called van Leer limiter, which 
actually has a functional behaviour quite close to our soft 
superbee. However, the basic numerical algorithm for solv- 
ing the hydrodynamic equation s differ consider a bly b etween 
our method and those used in iKlev fc Nelsoiil (|200'if ). so it 
is somewhat too simplified to compare them at the level of 
limiter functions. This clearly deserves more study. In the 
main text, we have studied both cases, using the minmod 
limiter to study a case for which the 3:1 Lindblad resonance 
does not operate. 



